Discovering Cyclic Causal Models

Supplementary Material

Author

Kyuri Park

Published

April 9, 2023


1 Introduction

Statistical network models have gained popularity in analyzing multivariate psychological data. These models often interpret network parameters as reflecting causal relationships, making it a form of causal discovery. However, recent research has shown that network models may not perform well as causal discovery tools for discovering acyclic causal structures (DAGs), and alternative methods are preferred for this task (Ryan, Bringmann, and Schuurman 2022).

But, acyclic causal models may not be suitable for some psychological phenomena, such as psychopathologies, which we expect to have cycles or feedback loop relationships between symptoms. While cyclic causal discovery methods have been developed in the computer science literature, they are not widely applied or well understood in empirical practice.

To address this gap, the main paper provides an accessible introduction to the basics of cyclic causal discovery for empirical researchers (Park 2023). It examines three different cyclic causal discovery methods and investigates their performance in typical psychological research contexts through a simulation study.

This supplementary material provides a more detailed analysis of the main simulation study. We omitted these results from the paper due to space constraints. In this supplementary material, we show the estimated PAGs and their associated measures of accuracy and uncertainty for each simulation condition. Additionally, we provide extra results on the running time of the algorithms used in the study.


2 Simulation Conditions

Below, we show the directed cyclic graph (DCG) and the corresponding true ancestral graph for each condition. There are in total eight conditions: model size (\(p = 5\), \(p = 10\)) \(\times\) density (sparse, dense) \(\times\) presence of latent confounder (presence, absence).

See code here.
#|label: simulation-specifics

## simulation design specifics
# specify the sample sizes
N <- c(50, 150, 500, 1000, 2000, 3000, 4000, 5000, 7500, 10000)
# specify replication number
n <- 500
# specify alpha level
alpha <- 0.01
# allow parallel processing
plan(multisession) 

2.1 5p_sparse

See code here.
#|label: 5psparse

## ====================
## 5p - sparse
## ====================

# specify B matrix
B5sparse = matrix(c(0, 0, 0, 0, 0,
                 1, 0, 0.8, 0, 0,
                 0, 0, 0, 0.9, 0,
                 0, 0.7, 0, 0, 1.5,
                 0, 0, 0, 0, 0), 5, 5, byrow = T)


colnames(B5sparse) <- c("X1", "X2", "X3", "X4", "X5")

# specify layout
layout5 = matrix(c(0,1,
                   0,0,
                   1,-1,
                   2,0,
                   2,1),5,2,byrow = T)

par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5psparse <- qgraph(t(B5sparse), layout=layout5, labels = colnames(B5sparse), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B5sparse)

# generate data 
simdata_5psparse <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B5sparse, N = z),  
            simplify = FALSE)
}, .options = furrr_options(seed=123))


## True Ancestral Graph
dcg_5psparse <- matrix(c(0,1,0,0,0,
                         0,0,0,1,0,
                         0,1,0,0,0,
                         0,0,1,0,0,
                         0,0,0,1,0), 5,5,byrow=T)
trueag_5psparse <- true_ancestral(dcg_5psparse, gen_dat(B5sparse), gaussCItest)
dimnames(trueag_5psparse) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5psparse)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
# Run CCD algorithm
# ccd_5psparse <- simdata_5psparse %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#       )
# mat_5psparse <- ccd_5psparse %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))

load("../data/mat_5psparse.RData")
# can we have a more elegant way to combine "graphNEL" plots? ? ? ? ? 
# par(mfrow=c(2,5))
# pag_ccd5psparse <- map2(ccd_5psparse, mat_5psparse,
#                         ~map2(.x, .y, plotPAG)
#                          ) 

## Run FCI algorithm
# fci_5psparse <- simdata_5psparse %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # extract amat
#       )
load("../data/fci_5psparse.RData")

## Run CCI algorithm
# cci_5psparse <- simdata_5psparse %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#       )
load("../data/cci_5psparse.RData")
# par(mfrow=c(2,5))
# pag_cci5psparse <- cci_5psparse %>%
#   map_depth(2, ~plotAG(.x))


## evaluation 
# CCD
res_ccd5psparse <- mat_5psparse %>% 
  map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd5psparse <- mat_5psparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% do.call("cbind", .) %>% apply(., 2, unlist) %>%  as.data.frame %>% 
  rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5psparse <- mat_5psparse %>% 
  map_depth(2, ~SHD(trueag_5psparse, .x)) %>% do.call("cbind", .) %>% apply(., 2, unlist) %>%  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci5psparse <- fci_5psparse %>% 
  map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci5psparse <- fci_5psparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5psparse <- fci_5psparse %>% 
  map_depth(2, ~SHD(trueag_5psparse, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci5psparse <- cci_5psparse %>% 
  map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci5psparse <- cci_5psparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5psparse <- cci_5psparse %>% 
  map_depth(2, ~SHD(trueag_5psparse, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.2 5p_dense

See code here.
#|label: 5pdense

## ====================
## 5p - dense
## ====================
# specify B matrix

## cyclic case
B5dense = matrix(c(0, 0, 0, 0, 0,
                    1.4, 0, 0.8, 0, 0,
                    0, 0, 0, 0.9, 0,
                    0, 0.7, 0, 0, 1,
                    1, 0, 0, 0, 0), 5, 5, byrow = T)


colnames(B5dense) <- c("X1", "X2", "X3", "X4", "X5")

par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5pdense <- qgraph(t(B5dense), layout=layout5, labels = colnames(B5dense), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B5dense)

# generate data (sample size as specified above)
simdata_5pdense <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B5dense, N = z),  
            simplify = FALSE)
}, .options = furrr_options(seed=123))


# True Ancestral Graph
dcg_5pdense <- matrix(c(0,1,0,0,1,
                        0,0,0,1,0,
                        0,1,0,0,0,
                        0,0,1,0,0,
                        0,0,0,1,0), 5,5,byrow=T)

trueag_5pdense <- true_ancestral(dcg_5pdense, gen_dat(B5dense), gaussCItest)
dimnames(trueag_5pdense) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5pdense)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_5pdense <- simdata_5pdense %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_5pdense <- ccd_5pdense %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))

# par(mfrow=c(2,5))
# pag_ccd5pdense <- map2(ccd_5pdense, mat_5pdense,
#                         ~map2(.x, .y, plotPAG)
#                          )
load("../data/mat_5pdense.RData")


## Run FCI algorithm
# fci_5pdense <- simdata_5pdense %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # extract amat
#   )
load("../data/fci_5pdense.RData")


## Run CCI algorithm
# cci_5pdense <- simdata_5pdense %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_5pdense.RData")

## evaluation
# CCD
res_ccd5pdense <- mat_5pdense %>% 
  map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame()
# UNCERTAINTY
uncer_ccd5pdense <- mat_5pdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pdense <- mat_5pdense %>% 
  map_depth(2, ~SHD(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci5pdense <- fci_5pdense %>% 
  map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci5pdense <- fci_5pdense%>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pdense <- fci_5pdense %>% 
  map_depth(2, ~SHD(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci5pdense <- cci_5pdense %>% 
  map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 

# UNCERTAINTY 
uncer_cci5pdense <- cci_5pdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pdense <- cci_5pdense %>% 
  map_depth(2, ~SHD(trueag_5pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.3 10p_sparse

See code here.
## ====================
## 10p - sparse
## ====================

B10sparse = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0.8, 0, 0, 0, 0, 0, 0, 0, 
                  0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                  0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 
                  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                  0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0, 
                  0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)

dimnames(B10sparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))

# specify layout
layout10 = matrix(c(0,1,
                      2,1,
                      1,0,
                      2,-1,
                      3,0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1),10,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10psparse <- qgraph(t(B10sparse), layout = layout10, labels = colnames(B10sparse), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B10sparse)

# generate data (sample size as specified above)
simdata_10psparse <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B10sparse, N = z),  
            simplify = FALSE)
}, .options = furrr_options(seed=123))

## True Ancestral Graph

dcg_10psparse <- matrix(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                          0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 
                          0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
                          0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
                          0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                          0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 
                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                          0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
                          0, 0, 0, 0, 0, 0, 0, 0, 1, 0), 10, 10, byrow = T)

trueag_10psparse <- true_ancestral(dcg_10psparse, gen_dat(B10sparse), gaussCItest)
dimnames(trueag_10psparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10psparse)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_10psparse <- simdata_10psparse %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# 
# mat_10psparse <- ccd_10psparse %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10psparse.RData")

## Run FCI algorithm
# fci_10psparse <- simdata_10psparse %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
load("../data/fci_10psparse.RData")

## Run CCI algorithm
# cci_10psparse <- simdata_10psparse %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_10psparse.RData")

## evaluation
# CCD
res_ccd10psparse <- mat_10psparse %>% 
  map_depth(2, 
  ~precision_recall(trueag_10psparse, .x)) %>%
  do.call("cbind", .) %>% t() %>%  
  apply(., 2, unlist) %>%  as.data.frame()
# UNCERTAINTY
uncer_ccd10psparse <- mat_10psparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10psparse <- mat_10psparse %>% 
  map_depth(2, ~SHD(trueag_10psparse, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci10psparse <- fci_10psparse %>% 
  map_depth(2, ~precision_recall(trueag_10psparse, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 

# UNCERTAINTY
uncer_fci10psparse <- fci_10psparse%>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10psparse <- fci_10psparse %>% 
  map_depth(2, ~SHD(trueag_10psparse, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci10psparse <- cci_10psparse %>% 
  map_depth(2, ~precision_recall(trueag_10psparse, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci10psparse <- cci_10psparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10psparse <- cci_10psparse %>% 
  map_depth(2, ~SHD(trueag_10psparse, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.4 10p_dense

See code here.
## ====================
## 10p - dense
## ====================

B10dense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     0.3, 0, 0.8, 0, 0, 0, 0, 0, 0, 0, 
                     0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                     0, 0, 1.1, 0, 0, 0.9, 0, 0, 0, 0, 
                     0, 0.4, 0, 1, 0, 0, 0, 0, 0, 0, 
                     0, 0, 0, 0, 0.9, 0, 0.5, 0, 0, 0, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0.8, 1, 
                     0, 0, 0, 0, 0, 0, 1, 0, 0, 0.4, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)



colnames(B10dense) <- c(paste("X", 1:10, sep=""))

# specify layout
layout10 = matrix(c(0,1,
                    2,1,
                    1,0,
                    2,-1,
                    3,0,
                    4, -1,
                    5, 0,
                    6, -1,
                    4, 1,
                    7, 1),10,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pdense <- qgraph(t(B10dense), layout = layout10, labels = colnames(B10dense), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B10dense)

# generate data (sample size as specified above)
simdata_10pdense <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B10dense, N = z),  
            simplify = FALSE)
}, .options = furrr_options(seed=123))


## True Ancestral Graph

dcg_10pdense <- matrix(c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
                         0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 
                         0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
                         0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
                         0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                         0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
                         0, 0, 0, 0, 0, 0, 1, 1, 1, 0), 10, 10, byrow = T)

trueag_10pdense <- true_ancestral(dcg_10pdense, gen_dat(B10dense), gaussCItest)
dimnames(trueag_10pdense) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10pdense)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_10pdense <- simdata_10pdense  %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_10pdense  <- ccd_10pdense  %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pdense.RData")


## Run FCI algorithm
# fci_10pdense <- simdata_10pdense  %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
load("../data/fci_10pdense.RData")

## Run CCI algorithm
# cci_10pdense  <- simdata_10pdense %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_10pdense.RData")


## evaluation
# CCD
res_ccd10pdense  <- mat_10pdense  %>% 
  map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd10pdense  <- mat_10pdense  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pdense  <- mat_10pdense %>% 
  map_depth(2, ~SHD(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))


# FCI
res_fci10pdense  <- fci_10pdense  %>% 
  map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci10pdense <- fci_10pdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pdense  <- fci_10pdense %>% 
  map_depth(2, ~SHD(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci10pdense <- cci_10pdense %>% 
  map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci10pdense <- cci_10pdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pdense <- cci_10pdense %>% 
  map_depth(2, ~SHD(trueag_10pdense, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.5 5p_LVsparse

See code here.
## ====================
## 5p with LV sparse
## ====================
# specify B matrix

B5_lvsparse = matrix(c(0, 0, 0, 0, 0, 1,
                 0, 0, 0.4, 0, 0, 1,
                 0, 0, 0, 0.5, 0,0,
                 0, 0.7, 0, 0, 1.5,0,
                 0, 0, 0, 0, 0,0,
                 0,0,0,0,0,0), 6, 6, byrow = T)

colnames(B5_lvsparse) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
                      0,0,
                      1,-1,
                      2,0,
                      2,1,
                      -1, 0.5),6,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5p_lvsparse <- qgraph(t(B5_lvsparse), layout=layout5_lv, labels = colnames(B5_lvsparse), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B5_lvsparse)

# generate data (sample size as specified above)
simdata_5pLVsparse <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B5_lvsparse, N = z)[,-6],  
            simplify = FALSE)
}, .options = furrr_options(seed=123))

## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_5psparseLV <- matrix(c(0, 2, 0, 0, 0, 
                           2, 0, 0, 1, 0, 
                           0, 1, 0, 0, 0,
                           0, 0, 1, 0, 0,
                           0, 0, 0, 1, 0), 5, 5, byrow = T)

trueag_5psparseLV <- true_ancestral(dcg_5psparseLV, gen_dat(B5_lvsparse), gaussCItest)
dimnames(trueag_5psparseLV) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5psparseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_5pLVsparse <- simdata_5pLVsparse  %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_5pLVsparse  <- ccd_5pLVsparse  %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_5pLVsparse.RData")

## Run FCI algorithm
# fci_5pLVsparse <- simdata_5pLVsparse  %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
load("../data/fci_5pLVsparse.RData")


## Run CCI algorithm
# cci_5pLVsparse  <- simdata_5pLVsparse %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_5pLVsparse.RData")



## evaluation
# CCD
res_ccd5pLVsparse  <- mat_5pLVsparse  %>% 
  map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd5pLVsparse  <- mat_5pLVsparse  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pLVsparse <- mat_5pLVsparse %>% 
  map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci5pLVsparse  <- fci_5pLVsparse  %>% 
  map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci5pLVsparse <- fci_5pLVsparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pLVsparse <- fci_5pLVsparse %>% 
  map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci5pLVsparse <- cci_5pLVsparse %>% 
  map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci5pLVsparse <- cci_5pLVsparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pLVsparse <- cci_5pLVsparse %>% 
  map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.6 5p_LVdense

See code here.
## ====================
## 5p with LV dense
## ====================
# specify B matrix

B5_lvdense = matrix(c(0, 0, 0, 0, 0, 1,
                       0, 0, 0.4, 0, 0, 1,
                       0, 0, 0, 0.5, 0,0,
                       0, 0.7, 0, 0, 0.5, 0,
                       0.6, 0, 0, 0, 0,0,
                       0,0,0,0,0,0), 6, 6, byrow = T)

colnames(B5_lvdense) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
                      0,0,
                      1,-1,
                      2,0,
                      2,1,
                      -1, 0.5),6,2,byrow = T)

par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5p_lvdense <- qgraph(t(B5_lvdense), layout=layout5_lv, labels = colnames(B5_lvdense), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)


## Data generating
# equilibrium check
equilibrium_check(B5_lvdense)

# generate data (sample size as specified above)
simdata_5pLVdense <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B5_lvdense, N = z)[,-6],  
            simplify = FALSE)
}, .options = furrr_options(seed=123))

## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_5pdenseLV <- matrix(c(0, 2, 0, 0, 1, 
                          2, 0, 0, 1, 0, 
                          0, 1, 0, 0, 0,
                          0, 0, 1, 0, 0,
                          0, 0, 0, 1, 0), 5, 5, byrow = T)

trueag_5pdenseLV <- true_ancestral(dcg_5pdenseLV, gen_dat(B5_lvdense), gaussCItest)
dimnames(trueag_5pdenseLV) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5pdenseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_5pLVdense <- simdata_5pLVdense  %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_5pLVdense  <- ccd_5pLVdense  %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_5pLVdense.RData")


## Run FCI algorithm
# fci_5pLVdense <- simdata_5pLVdense  %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#                     alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
# 
load("../data/fci_5pLVdense.RData")

## Run CCI algorithm
# cci_5pLVdense  <- simdata_5pLVdense %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
# 
load("../data/cci_5pLVdense.RData")


## evaluation
# CCD
res_ccd5pLVdense  <- mat_5pLVdense  %>% 
  map_depth(2, ~precision_recall(trueag_5pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd5pLVdense  <- mat_5pLVdense  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pLVdense <- mat_5pLVdense %>% 
  map_depth(2, ~SHD(trueag_5pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci5pLVdense  <- fci_5pLVdense  %>% 
  map_depth(2, ~precision_recall(trueag_5pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci5pLVdense <- fci_5pLVdense  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pLVdense  <- fci_5pLVdense  %>% 
  map_depth(2, ~SHD(trueag_5pdenseLV , .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci5pLVdense  <- cci_5pLVdense  %>% 
  map_depth(2, ~precision_recall(trueag_5pdenseLV , .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci5pLVdense  <- cci_5pLVdense  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pLVdense  <- cci_5pLVdense  %>% 
  map_depth(2, ~SHD(trueag_5pdenseLV , .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.7 10p_LVsparse

See code here.
## ====================
## 10p with LV sparse
## ====================
# specify B matrix

## with 2 LVs
B10_lvsparse = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 0.7,
                        0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.7,
                        0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.8, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 12, 12, byrow = T)
colnames(B10_lvsparse) <- c(paste("X", 1:10, sep=""), "L1", "L2")


# specify layout
layout10LV2 = matrix(c(0, 1,
                      2, 1,
                      1, 0,
                      2, -1,
                      3, 0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1,
                      8, 0,
                      3, 2), 12, 2, byrow = T)

par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pLVsparse <- qgraph(t(B10_lvsparse), layout = layout10LV2, labels = colnames(B10_lvsparse), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B10_lvsparse)

# generate data (sample size as specified above)
simdata_10pLVsparse <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B10_lvsparse, N = z)[,-c(11,12)],  
            simplify = FALSE)
}, .options = furrr_options(seed=123))

## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_10psparseLV <- matrix(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 
                            0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 
                            0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
                            0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 
                            0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                            0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 
                            0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 
                            0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
                            0, 0, 0, 0, 0, 0, 0, 2, 1, 0), 10, 10, byrow = T)

trueag_10psparseLV <- true_ancestral(dcg_10psparseLV, gen_dat(B10_lvsparse), gaussCItest)
dimnames(trueag_10psparseLV) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10psparseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_10pLVsparse  <- simdata_10pLVsparse   %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_10pLVsparse   <- ccd_10pLVsparse %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pLVsparse.RData")


## Run FCI algorithm
# fci_10pLVsparse  <- simdata_10pLVsparse   %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
load("../data/fci_10pLVsparse.RData")

## Run CCI algorithm
# cci_10pLVsparse  <- simdata_10pLVsparse %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_10pLVsparse.RData")


## evaluation
# CCD
res_ccd10pLVsparse   <- mat_10pLVsparse  %>% 
  map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd10pLVsparse  <- mat_10pLVsparse  %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pLVsparse <- mat_10pLVsparse %>% 
  map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# FCI
res_fci10pLVsparse <- fci_10pLVsparse  %>% 
  map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci10pLVsparse <- fci_10pLVsparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pLVsparse <- fci_10pLVsparse %>% 
  map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI 
res_cci10pLVsparse <- cci_10pLVsparse %>% 
  map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci10pLVsparse <- cci_10pLVsparse %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pLVsparse <- cci_10pLVsparse %>% 
  map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

2.8 10p_LVdense

See code here.
## ====================
## 10p with LV dense
## ====================
# specify B matrix
B10_lvdense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0.5, 0, 1.4, 0, 0, 0, 0, 0, 0, 0, 0, 1.1,
                       0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 1, 0, 0, 0.9, 0, 0, 0, 0, 0, 0,
                       0, 0, 1.2, 0.8, 0, 0, 0, 0, 0, 0, 0, 0.7,
                       0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 1.4, 0, 0.4, 0,
                       0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.6, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 12, 12, byrow = T)

colnames(B10_lvdense) <- c(paste("X", 1:10, sep=""), "L1", "L2")

# specify layout
layout10LV2 = matrix(c(0, 1,
                      2, 1,
                      1, 0,
                      2, -1,
                      3, 0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1,
                      8, 0,
                      3, 2), 12, 2, byrow = T)

par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pLVdense <- qgraph(t(B10_lvdense), layout = layout10LV2, labels = colnames(B10_lvdense), theme="colorblind")
title("True cyclic graph",  font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)

## Data generating
# equilibrium check
equilibrium_check(B10_lvdense)

# generate data (sample size as specified above)
simdata_10pLVdense <- N %>% future_map(function(z) {
  replicate(n = n,
            expr = gen_dat(B10_lvdense, N = z)[,-c(11,12)],  
            simplify = FALSE)
}, .options = furrr_options(seed=123))


## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_10pdenseLV <- matrix(c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
                           0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 
                           0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 
                           0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
                           0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 
                           0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                           0, 0, 0, 0, 0, 1, 0, 1, 0, 2, 
                           0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 
                           0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
                           0, 0, 0, 0, 0, 0, 2, 2, 1, 0), 10, 10, byrow = T)


trueag_10pdenseLV <- true_ancestral(dcg_10pdenseLV, gen_dat(B10_lvdense), gaussCItest)
dimnames(trueag_10pdenseLV) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10pdenseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)

See code here.
## Run CCD algorithm
# ccd_10pLVdense  <- simdata_10pLVdense   %>%
#   map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
#   )
# mat_10pLVdense   <- ccd_10pLVdense %>% 
#   map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pLVdense.RData")

## Run FCI algorithm
# fci_10pLVdense  <- simdata_10pLVdense   %>%
#   map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
#            alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
#   )
load("../data/fci_10pLVdense.RData")

## Run CCI algorithm
# cci_10pLVdense  <- simdata_10pLVdense  %>%
#   map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag  # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
#   )
load("../data/cci_10pLVdense.RData")


## evaluation
# CCD
res_ccd10pLVdense   <- mat_10pLVdense  %>% 
  map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_ccd10pLVdense  <- mat_10pLVdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pLVdense <- mat_10pLVdense %>% 
  map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N)) 

# FCI
res_fci10pLVdense <- fci_10pLVdense  %>% 
  map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_fci10pLVdense <- fci_10pLVdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pLVdense <- fci_10pLVdense %>% 
  map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

# CCI
res_cci10pLVdense <- cci_10pLVdense %>% 
  map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% t() %>%  apply(., 2, unlist) %>%  as.data.frame() 
# UNCERTAINTY
uncer_cci10pLVdense <- cci_10pLVdense %>% 
  map_depth(2, ~uncertainty(.x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pLVdense <- cci_10pLVdense %>% 
  map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>% 
  do.call("cbind", .) %>% apply(., 2, unlist) %>%  
  as.data.frame %>% rename_with(~ paste0("N = ", N))

3 Performance Evaluation per Condition

Below, we present the most frequently estimated PAG for each sample size (N) from each algorithm. The graphs are obtained by picking the most frequent type of edge-endpoint produced by algorithms out of 500 simulations. Additionally, we present the matrix plots of the correct estimation and uncertainty for each condition. The darker the color (blue/red) is, the higher the rate of correct estimation (blue) or uncertainty (red).

See code here.
par(mfrow=c(2,5))
vec <- list("CCD-5p-sparse" = mat_5psparse, "FCI-5p-sparse" = fci_5psparse, "CCI-5p-sparse" = cci_5psparse)

graphs_5psparse <- list()
for(i in seq_along(vec)){
  ## high-freq graphs
  graphs_5psparse[[i]] <- vec[[i]] %>% 
  map(~high_freq(.x, p = 5) %>% 
        plotAG)
  ## correct prop plots
  vec[[i]] %>% 
    imap(
      ~ prop_correct(.x, trueag_5psparse, p = 5) %>% 
        # long format
        reshape2::melt() %>% 
        ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
        geom_tile() +
        geom_text() + 
        # reverse factor level
        scale_y_discrete(limits=rev) + 
        scale_fill_gradient(low="grey90", high="blue") +
        labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
    ) %>% 
    ggpubr::ggarrange(plotlist = .,
                      ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
    ggpubr::annotate_figure(top = text_grob(glue::glue("Correct Proportion {names(vec[[i]])}"), face = "bold", size = 18, family = "Palatino"))
  ## uncertainty prop plots
  vec[[i]] %>% 
    imap(
      ~ prop_uncertain(.x, p = 5) %>% 
        # long format
        reshape2::melt() %>% 
        ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
        geom_tile() +
        geom_text() + 
        # reverse factor level
        scale_y_discrete(limits=rev) + 
        scale_fill_gradient(low="grey90", high="red") +
        labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
    ) %>% 
    ggpubr::ggarrange(plotlist = .,
                      ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
    ggpubr::annotate_figure(top = text_grob(glue::glue("Uncertainty {names(vec)[i]}"), face = "bold", size = 18, family = "Palatino"))
}

corprop_plots_5psparse <- list()
# prop correct
for(i in seq_along(vec)){
  corprop_plots_5psparse[[i]] <- vec[[i]] %>% 
    imap(
      ~ prop_correct(.x, trueag_5psparse, p = 5) %>% 
        # long format
        reshape2::melt() %>% 
        ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
        geom_tile() +
        geom_text() + 
        # reverse factor level
        scale_y_discrete(limits=rev) + 
        scale_fill_gradient(low="grey90", high="blue") +
        labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
    ) %>% 
    ggpubr::ggarrange(plotlist = .,
                      ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
    ggpubr::annotate_figure(top = text_grob(glue::glue("Correct Proportion {names(vec[[i]])}"), face = "bold", size = 18, family = "Palatino"))
}

ucprop_plots_5psparse <- list()
# prop uncertain
for(i in seq_along(vec)){
  ucprop_plots_5psparse[[i]] <- vec[[i]] %>% 
    imap(
      ~ prop_uncertain(.x, p = 5) %>% 
        # long format
        reshape2::melt() %>% 
        ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
        geom_tile() +
        geom_text() + 
        # reverse factor level
        scale_y_discrete(limits=rev) + 
        scale_fill_gradient(low="grey90", high="red") +
        labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
    ) %>% 
    ggpubr::ggarrange(plotlist = .,
                      ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
    ggpubr::annotate_figure(top = text_grob(glue::glue("Uncertainty {names(vec)[i]}"), face = "bold", size = 18, family = "Palatino"))
}
See code here.
## CCD 5p sparse case
# high frequency
par(mfrow=c(2,5))
mat_5psparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
mat_5psparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparse, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_5psparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## FCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

fci_5psparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
fci_5psparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparse, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_5psparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

cci_5psparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
cci_5psparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparse, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_5psparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 5p dense case
# high frequency 
par(mfrow=c(2,5))
mat_5pdense %>% 
  map(~high_freq(.x, p = 5) %>% 
        plotAG)

See code here.
# prop correct
mat_5pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_5pdense, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_5pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## FCI 5p dense case
# high frequency 
par(mfrow=c(2,5))

fci_5pdense %>% 
  map(~high_freq(.x, p = 5) %>% 
        plotAG)

See code here.
# prop correct
fci_5pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_5pdense, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_5pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p dense case
# high frequency 
par(mfrow=c(2,5))

cci_5pdense %>% 
  map(~high_freq(.x, p = 5) %>% 
        plotAG)

See code here.
# prop correct
cci_5pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_5pdense, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_5pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 5) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 10p sparse case
# high frequency 
par(mfrow=c(2,5))
mat_10psparse %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
mat_10psparse %>% 
  imap(
    ~ prop_correct(.x, trueag_10psparse, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_10psparse %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## FCI 10p sparse case
# high frequency 
par(mfrow=c(2,5))

fci_10psparse %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
fci_10psparse %>% 
  imap(
    ~ prop_correct(.x, trueag_10psparse, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_10psparse %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 10p sparse case
# high frequency 
par(mfrow=c(2,5))

cci_10psparse %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
cci_10psparse %>% 
  imap(
    ~ prop_correct(.x, trueag_10psparse, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_10psparse %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 10p dense case
# high frequency 
par(mfrow=c(2,5))
mat_10pdense %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
mat_10pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_10pdense, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_10pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## FCI 10p dense case
# high frequency 
par(mfrow=c(2,5))

fci_10pdense %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
fci_10pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_10psparse, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_10pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 10p dense case
# high frequency 
par(mfrow=c(2,5))

cci_10pdense %>% 
  map(~high_freq(.x, p = 10) %>% 
        plotAG)

See code here.
# prop correct
cci_10pdense %>% 
  imap(
    ~ prop_correct(.x, trueag_10psparse, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_10pdense %>% 
  imap(
    ~ prop_uncertain(.x, p = 10) %>% 
      # long format
      reshape2::melt() %>% 
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) + 
      geom_tile() +
      geom_text() + 
      # reverse factor level
      scale_y_discrete(limits=rev) + 
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}")) 
  ) %>% 
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 5p sparse LV case
# high frequency
par(mfrow=c(2,5))
mat_5pLVsparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
mat_5pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_5pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# high frequency
par(mfrow=c(2,5))

fci_5pLVsparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
fci_5pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_5pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

cci_5pLVsparse %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
cci_5pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_5pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 5p dense LV case
# high frequency
par(mfrow=c(2,5))
mat_5pLVdense %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
mat_5pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_5pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# high frequency
par(mfrow=c(2,5))

fci_5pLVdense %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
fci_5pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_5pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

cci_5pLVdense %>%
  map(~high_freq(.x, p = 5) %>%
        plotAG)

See code here.
# prop correct
cci_5pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_5pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 5) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 10p sparse LV case
# high frequency
par(mfrow=c(2,5))
mat_10pLVsparse %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
mat_10pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_10pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# high frequency
par(mfrow=c(2,5))

fci_10pLVsparse %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
fci_10pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_10pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

cci_10pLVsparse %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
cci_10pLVsparse %>%
  imap(
    ~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_10pLVsparse %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCD 10p dense LV case
# high frequency
par(mfrow=c(2,5))
mat_10pLVdense %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
mat_10pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
mat_10pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# high frequency
par(mfrow=c(2,5))

fci_10pLVdense %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
fci_10pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
fci_10pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))

cci_10pLVdense %>%
  map(~high_freq(.x, p = 10) %>%
        plotAG)

See code here.
# prop correct
cci_10pLVdense %>%
  imap(
    ~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="blue") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))

See code here.
# prop uncertain
cci_10pLVdense %>%
  imap(
    ~ prop_uncertain(.x, p = 10) %>%
      # long format
      reshape2::melt() %>%
      ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
      geom_tile() +
      geom_text() +
      # reverse factor level
      scale_y_discrete(limits=rev) +
      scale_fill_gradient(low="grey90", high="red") +
      labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
  ) %>%
  ggpubr::ggarrange(plotlist = .,
                    ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
  ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))


4 Overall Performance Evalulation

Here we summarize the performance of the algorithms across all conditions using the evaluation metrics: structural Hamming distance (SHD), precision, recall, and uncertainty. Each point represents the average value of each metric and shaded area represents interquartile range (IQR).

See code here.
## ============================
## Create neat dataframe
## ============================

## Compute average precision & recall and corresponding sd for each condition
pre_rec <- list(
  # put all the results together in a list
  res_ccd5psparse, res_fci5psparse, res_cci5psparse, res_ccd10psparse, res_fci10psparse, res_cci10psparse, res_ccd5pdense, res_fci5pdense, res_cci5pdense, res_ccd10pdense, res_fci10pdense, res_cci10pdense, res_ccd5pLVsparse, res_fci5pLVsparse, res_cci5pLVsparse, res_ccd5pLVdense, res_fci5pLVdense, res_cci5pLVdense, res_ccd10pLVsparse, res_fci10pLVsparse, res_cci10pLVsparse, res_ccd10pdense,  res_fci10pLVdense, res_cci10pLVdense
) %>% 
  # transpose df
  map(~ sjmisc::rotate_df(.x) %>%
        # add sample size (N) info
        rename_with(~paste0(.x, "N = ", rep(N, each=8)))  %>%
        # think about how to deal with NAs or do I want to define sth. else instead of NAs.
        # na.omit(.x) %>% 
        # get the average and sd
        dplyr::summarise(across(everything(.), list(mean = ~mean(., na.rm=T), sd = ~sd(., na.rm=T))))) %>% 
  bind_rows() %>% 
  mutate(algorithm = rep(c("CCD", "FCI", "CCI"), 8),
         condition = rep(c("5p_sparse", "10p_sparse", "5p_dense", "10p_dense", "5p_LVsparse", "5p_LVdense", "10p_LVsparse", "10p_LVdense"), each=3),
         netsize = stringr::str_split(condition, "_", simplify=T)[,1],
         latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
         densities = stringr::str_remove(stringr::str_split(condition, "_", simplify=T)[,2], "LV")
  ) %>%
  # brings the algorithm and condition names first
  relocate(where(is.character), .before = where(is.numeric)) %>% 
  # convert it to a long format
  tidyr::pivot_longer(!c(algorithm, condition, netsize, latentvar, densities), names_to = "metric", values_to = "value") %>% 
  # Add sample size column (N) & clean up the column name 
  mutate(N = stringr::str_extract(metric, "(?<=[N =])\\d+"),
         metric = stringr::str_replace_all(metric, "[0-9.]+|[N =]", "")) 


## Compute average uncertainty rate and corresponding sd for each condition
uncertainties <- bind_rows(
  # bind all results from each condition
  "CCD_5p-sparse" = uncer_ccd5psparse, "FCI_5p-sparse" = uncer_fci5psparse, "CCI_5p-sparse"=uncer_cci5psparse, "CCD_10p-sparse"=uncer_ccd10psparse, "FCI_10p-sparse" = uncer_fci10psparse, "CCI_10p-sparse" = uncer_cci10psparse, "CCD_5p-dense"=uncer_ccd5pdense, "FCI_5p-dense"=uncer_fci5pdense, "CCI_5p-dense"=uncer_cci5pdense, "CCD_10p-dense"=uncer_ccd10pdense, "FCI_10p-dense"=uncer_fci10pdense, "CCI_10p-dense"=uncer_cci10pdense, "CCD_5p-LVsparse"=uncer_ccd5pLVsparse, "FCI_5p-LVsparse"=uncer_fci5pLVsparse, "CCI_5p-LVsparse"=uncer_cci5pLVsparse, "CCD_10p-LVsparse"=uncer_ccd10pLVsparse, "FCI_10p-LVsparse"=uncer_fci10pLVsparse, "CCI_10p-LVsparse"=uncer_cci10pLVsparse,
  "CCD_5p-LVdense"=uncer_ccd5pLVdense, "FCI_5p-LVdense"=uncer_fci5pLVdense, "CCI_5p-LVdense"=uncer_cci5pLVdense, "CCD_10p-LVdense"=uncer_ccd10pLVdense, "FCI_10p-LVdense"=uncer_fci10pLVdense, "CCI_10p-LVdense"=uncer_cci10pLVdense, .id="id"
) %>% 
  group_by(id) %>% 
  # get the average and sd
  summarise_all(list(means = mean, sds = sd)) %>%  
  mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
         condition = stringr::str_split(id, "_", simplify = T)[,2],
         netsize = stringr::str_split(condition, "-", simplify=T)[,1],
         latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
         densities = stringr::str_remove(stringr::str_split(condition, "-", simplify=T)[,2], "LV")
  ) %>% 
  tidyr::pivot_longer(!c(algorithm, condition, id, netsize, latentvar, densities), names_to = "name", values_to = "value") %>% 
  mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
         statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>% 
  dplyr::select(-id, -name) %>%  relocate(where(is.character), .before = where(is.numeric))


## Compute average SHD values and corresponding sd for each condition
SHDs <- bind_rows(
  # bind all results from each condition
  "CCD_5p-sparse" = SHD_ccd5psparse, "FCI_5p-sparse" = SHD_fci5psparse, "CCI_5p-sparse"=SHD_cci5psparse, "CCD_10p-sparse"= SHD_ccd10psparse, "FCI_10p-sparse" = SHD_fci10psparse, "CCI_10p-sparse" = SHD_cci10psparse, "CCD_5p-dense"= SHD_ccd5pdense, "FCI_5p-dense"=SHD_fci5pdense, "CCI_5p-dense"=SHD_cci5pdense, "CCD_10p-dense"= SHD_ccd10pdense, "FCI_10p-dense"=SHD_fci10pdense, "CCI_10p-dense"=SHD_cci10pdense, "CCD_5p-LVsparse"=SHD_ccd5pLVsparse, "FCI_5p-LVsparse"=SHD_fci5pLVsparse, "CCI_5p-LVsparse"=SHD_cci5pLVsparse, "CCD_10p-LVsparse"=SHD_ccd10pLVsparse, "FCI_10p-LVsparse"=SHD_fci10pLVsparse, "CCI_10p-LVsparse"=SHD_cci10pLVsparse, 
  "CCD_5p-LVdense"=SHD_ccd5pLVdense, "FCI_5p-LVdense"=SHD_fci5pLVdense, "CCI_5p-LVdense"=SHD_cci5pLVdense, "CCD_10p-LVdense"=SHD_ccd10pLVdense, "FCI_10p-LVdense"=SHD_fci10pLVdense, "CCI_10p-LVdense"=SHD_cci10pLVdense, .id="id"
) %>% 
  group_by(id) %>% 
  # get the average and sd
  summarise_all(list(means = mean, sds = sd)) %>%  
  mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
         condition = stringr::str_split(id, "_", simplify = T)[,2],
         netsize = stringr::str_split(condition, "-", simplify=T)[,1],
         latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
         densities = stringr::str_remove(stringr::str_split(condition, "-", simplify=T)[,2], "LV")
  ) %>% 
  tidyr::pivot_longer(!c(algorithm, condition, id, netsize, latentvar, densities), names_to = "name", values_to = "value") %>% 
  mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
         statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>% 
  dplyr::select(-id, -name) %>%  relocate(where(is.character), .before = where(is.numeric)) 




## ============================
## Create figures
## ============================

## specify the common figure theme
MyTheme <-  theme(plot.title = element_text(face = "bold", family = "Palatino", size = 15, hjust=0.5),
                  plot.subtitle = element_text(face = "italic", family = "Palatino", size = 15, hjust=0.5),
                  axis.text=element_text(face = "bold",family = "Palatino", size = 11),
                  axis.text.x = element_text(angle = 45, hjust = 1.2, vjust =1.2),
                  axis.title = element_text(face = "bold",family = "Palatino", size = 12),
                  legend.text = element_text(face = "bold", family = "Palatino", size = 12),
                  legend.position="bottom",
                  strip.text = element_text(face="bold", size=13, family = "Palatino"),
                  strip.background = element_rect(fill="#f0f0f0", linetype = "solid", color="gray"),
                  strip.placement = "outside",
                  panel.border = element_rect(color = "#DCDCDC", fill = NA),
                  panel.spacing.y = unit(4, "mm")
)


## SHD figure
SHDs %>%
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= as.numeric(N), y=means, group = algorithm, colour = algorithm, fill = algorithm)) +
  # add line plots
  geom_line(aes(group = algorithm)) +
  # add scattered points
  geom_point(size=1) + 
  # add interquartile range (IQR)
  geom_ribbon(aes(ymin=means+qnorm(0.25)*sds, ymax=means+qnorm(0.75)*sds), alpha=0.15, color=NA) +
  # specify custom colors
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "") +
  # apply the theme
  theme_minimal() +
  MyTheme + 
  # create a facet
  ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels = c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")),  scales = "free_y", switch="y") +
  scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
  ggtitle("SHD") 

See code here.
## Precision figure
pre_rec %>% 
  filter(grepl("average_precision", metric)) %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= as.numeric(N), y=average_precision_mean, group = algorithm, colour = algorithm, fill=algorithm)) +
  # add line plots
  geom_line(aes(group = algorithm)) +
  # add scattered points
  geom_point(size=1) +
  # add interquartile range (IQR)
  geom_ribbon(aes(ymin=average_precision_mean+qnorm(0.25)*average_precision_sd, ymax=average_precision_mean+qnorm(0.75)*average_precision_sd), alpha=0.15, color=NA) +
  # specify custom colors
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  # apply the theme
  theme_minimal() +
  MyTheme + 
  # create a facet
  ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")),  switch="y") +
  scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
  labs(title = "Precision", x = "", y = "")

See code here.
## Recall figure
pre_rec %>% 
  filter(grepl("average_recall", metric)) %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= as.numeric(N), y=average_recall_mean, group = algorithm, colour = algorithm, fill= algorithm)) +
  # add line plots
  geom_line(aes(group = algorithm)) +
  # add scattered points
  geom_point(size=1) +
  # add interquartile range (IQR)
  geom_ribbon(aes(ymin=average_recall_mean+qnorm(0.25)*average_recall_sd, ymax=average_recall_mean+qnorm(0.75)*average_recall_sd), alpha=0.15, color=NA) +
  # specify custom colors
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  # apply the theme
  theme_minimal() +
  MyTheme + 
  # create a facet
  ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")),  switch="y") +
  scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
  labs(title = "Recall", x = "", y = "")

See code here.
## Uncertainty figure
uncertainties %>%
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= as.numeric(N), y=means, group = algorithm, colour = algorithm, fill = algorithm)) +
  # add line plots
  geom_line(aes(group = algorithm)) +
  # add scattered points
  geom_point(size=1) + 
  # add interquartile range (IQR)
  geom_ribbon(aes(ymin=means+qnorm(0.25)*sds, ymax=means+qnorm(0.75)*sds), alpha=0.15, color=NA) +
  # specify custom colors
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "") +
  # apply the theme
  theme_minimal() +
  MyTheme + 
  # create a facet
  ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")),  scales = "free_y", switch="y") +
  scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
  ggtitle("Uncertainty")

5 Extra: Algorithm Running Time

The running time of the algorithms, presented in log(ms), indicates that FCI generally takes the shortest time, while CCD takes the longest across all conditions.

See code here.
# specify sig. level
alpha <- 0.05

times <- microbenchmark(
  ccd_5psparse = ccdKP(df=simdata_5psparse[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_5psparse = fci(list(C = cor(simdata_5psparse[[10]][[1]]), n = 1e1), indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5psparse[[10]][[1]])),
  cci_5psparse = cci(list(C = cor(simdata_5psparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5psparse[[10]][[1]])),

  ccd_5pdense = ccdKP(df=simdata_5pdense[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_5pdense = fci(list(C = cor(simdata_5pdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pdense[[10]][[1]])),
  cci_5pdense = cci(list(C = cor(simdata_5pdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pdense[[10]][[1]])),

  ccd_10psparse = ccdKP(df=simdata_10psparse[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_10psparse = fci(list(C = cor(simdata_10psparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(simdata_10psparse[[10]][[1]])),
  cci_10psparse = cci(list(C = cor(simdata_10psparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10psparse[[10]][[1]])),

  ccd_10pdense = ccdKP(df=simdata_10pdense[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_10pdense = fci(list(C = cor(simdata_10pdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pdense[[10]][[1]])),
  cci_10pdense = cci(list(C = cor(simdata_10pdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pdense[[10]][[1]])),

  ccd_5pLVsparse = ccdKP(df=simdata_5pLVsparse[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_5pLVsparse = fci(list(C = cor(simdata_5pLVsparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
                alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pLVsparse[[10]][[1]])),
  cci_5pLVsparse = cci(list(C = cor(simdata_5pLVsparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pLVsparse[[10]][[1]])),

  ccd_5pLVdense = ccdKP(df=simdata_5pLVdense[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_5pLVdense = fci(list(C = cor(simdata_5pLVdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
                       alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pLVdense[[10]][[1]])),
  cci_5pLVdense = cci(list(C = cor(simdata_5pLVdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pLVdense[[10]][[1]])),
  
  ccd_10pLVsparse = ccdKP(df=simdata_10pLVsparse[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_10pLVsparse = fci(list(C = cor(simdata_10pLVsparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
                alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pLVsparse[[10]][[1]])),
  cci_10pLVsparse = cci(list(C = cor(simdata_10pLVsparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pLVsparse[[10]][[1]])),
  
  ccd_10pLVdense = ccdKP(df=simdata_10pLVdense[[10]][[1]], dataType = "continuous", alpha = alpha),
  fci_10pLVdense = fci(list(C = cor(simdata_10pLVdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
                        alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pLVdense[[10]][[1]])),
  cci_10pLVdense = cci(list(C = cor(simdata_10pLVdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pLVdense[[10]][[1]]))
)

times <- times %>%
  mutate(algorithm = substr(expr, 1, 3),
         condition = stringr::str_split(expr, "_", simplify=T)[,2])

## plot the results
times %>%
  ggplot(aes(x=factor(condition, levels= c("5psparse", "5pdense", "10psparse", "10pdense", "5pLVsparse","5pLVdense", "10pLVsparse","10pLVdense")), y = log(time), col= factor(algorithm))) +
  geom_boxplot(position = "dodge",   outlier.size = 0.8, outlier.alpha = 0.2) + theme_classic() +
  # scale_x_discrete(name ="Condition",
  #                  labels=c("", "5p-sparse", "", "","5p-dense","","", "10p-sparse","","","10p-dense","","","5p-LV","","","10p-LV","")) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(y = " log(ms)", x = "", title = "Algorithm Running Time", subtitle = "Time in milliseconds (ms)") +
  theme(axis.text.x = element_text(face = "bold", family =  "Palatino", margin = margin(t = 13), size=10),
        legend.position="bottom",
        plot.subtitle = element_text(face = "italic", family = "Palatino"),
        plot.title = element_text(family = "Palatino", size=15),
        axis.title = element_text(family = "Palatino"),
        legend.text = element_text(family =  "Palatino"))


6 References

Park, Kyuri. 2023. Discovering cyclic causal models in psychological research.” GitHub repository. https://github.com/KyuriP/Thesis_KP.
Ryan, Oisín, Laura F. Bringmann, and Noémi K. Schuurman. 2022. “The Challenge of Generating Causal Hypotheses Using Network Models.” Structural Equation Modeling: A Multidisciplinary Journal 0 (0): 1–18. https://doi.org/10.1080/10705511.2022.2056039.

7 Session Information

sessionInfo()
R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.0

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] microbenchmark_1.4.9 CCI.KP_0.1.0         MASS_7.3-58.3       
 [4] ggpubr_0.6.0         ggplot2_3.4.2        qgraph_1.9.4        
 [7] dplyr_1.1.1          furrr_0.3.1          future_1.32.0       
[10] purrr_1.0.1          magrittr_2.0.3       Rgraphviz_2.40.0    
[13] graph_1.74.0         BiocGenerics_0.42.0  pcalg_2.7-8         
[16] DOT_0.1              rcausal_1.2.1        devtools_2.4.5      
[19] usethis_2.1.6        rJava_1.0-6         

loaded via a namespace (and not attached):
  [1] backports_1.4.1   Hmisc_5.0-1       plyr_1.8.8        igraph_1.4.1     
  [5] listenv_0.9.0     fastICA_1.2-3     digest_0.6.31     htmltools_0.5.5  
  [9] fansi_1.0.4       checkmate_2.1.0   memoise_2.0.1     cluster_2.1.4    
 [13] sfsmisc_1.1-14    remotes_2.4.2     globals_0.16.2    bdsmatrix_1.3-6  
 [17] prettyunits_1.1.1 ggh4x_0.2.4       jpeg_0.1-10       colorspace_2.1-0 
 [21] xfun_0.38         callr_3.7.3       crayon_1.5.2      jsonlite_1.8.4   
 [25] glue_1.6.2        gtable_0.3.3      V8_4.2.2          sjmisc_2.8.9     
 [29] car_3.1-2         pkgbuild_1.4.0    DEoptimR_1.0-11   ggm_2.5          
 [33] abind_1.4-5       scales_1.2.1      rstatix_0.7.2     miniUI_0.1.1.1   
 [37] Rcpp_1.0.10       xtable_1.8-4      htmlTable_2.4.1   clue_0.3-64      
 [41] foreign_0.8-84    Formula_1.2-5     stats4_4.2.3      profvis_0.3.7    
 [45] htmlwidgets_1.6.2 lavaan_0.6-15     ellipsis_0.3.2    urlchecker_1.0.1 
 [49] pkgconfig_2.0.3   farver_2.1.1      nnet_7.3-18       utf8_1.2.3       
 [53] tidyselect_1.2.0  labeling_0.4.2    rlang_1.1.0       reshape2_1.4.4   
 [57] later_1.3.0       munsell_0.5.0     tools_4.2.3       cachem_1.0.7     
 [61] cli_3.6.1         generics_0.1.3    sjlabelled_1.2.0  broom_1.0.4      
 [65] fdrtool_1.2.17    evaluate_0.20     stringr_1.5.0     fastmap_1.1.1    
 [69] yaml_2.3.7        processx_3.8.0    knitr_1.42        fs_1.6.1         
 [73] robustbase_0.95-1 glasso_1.11       pbapply_1.7-0     RBGL_1.72.0      
 [77] nlme_3.1-162      mime_0.12         compiler_4.2.3    rstudioapi_0.14  
 [81] curl_5.0.0        png_0.1-8         ggsignif_0.6.4    tibble_3.2.1     
 [85] pbivnorm_0.6.0    stringi_1.7.12    ps_1.7.4          lattice_0.21-8   
 [89] Matrix_1.5-4      psych_2.3.3       vctrs_0.6.1       pillar_1.9.0     
 [93] lifecycle_1.0.3   data.table_1.14.8 cowplot_1.1.1     insight_0.19.1   
 [97] corpcor_1.6.10    httpuv_1.6.9      R6_2.5.1          promises_1.2.0.1 
[101] gridExtra_2.3     parallelly_1.35.0 sessioninfo_1.2.2 codetools_0.2-19 
[105] gtools_3.9.4      pkgload_1.3.2     withr_2.5.0       mnormt_2.1.1     
[109] parallel_4.2.3    quadprog_1.5-8    rpart_4.1.19      tidyr_1.3.0      
[113] rmarkdown_2.21    carData_3.0-5     shiny_1.7.4       base64enc_0.1-3